Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données et produire un rapport.
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Notre jeu de données comporte 102 individus décrits par 6 variables :
head(data) %>% paged_table()
viz_edu <- data %>%
ggplot(aes(x = prestige, y = education, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(6, 16, 2)
)
viz_edu <- ggplotly(viz_edu)
viz_inc <- data %>%
ggplot(aes(x = prestige, y = income, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 20000, 10000)
)
viz_inc <- ggplotly(viz_inc)
viz_women <- data %>%
ggplot(aes(x = prestige, y = women, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 100, 25)
)
viz_women <- ggplotly(viz_women)
viz_census <- data %>%
ggplot(aes(x = prestige, y = census, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(2500, 7500, 2500)
)
viz_census <- ggplotly(viz_census)
subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)
En statistiques et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.
Dans un premier temps, nous pouvons étudier rapidement d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.
generate_scatterM <- function(data, mapping){
data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) +
geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}
scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))
scatterM <- ggplotly(scatterM)
scatterM
Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot pour une meilleure représentation.
education <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(20, 80, 20)
)
scatter_edu <- ggplotly(education)
scatter_edu
La droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.
La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.
On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.
Nous testons aussi la relation entre les variables income et prestige.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(0, 25000, 5000)
) +
scale_y_continuous(
breaks = seq(20, 100, 20)
)
scatter_inc <- ggplotly(income)
scatter_inc
Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.
Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.
lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu
##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361
« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.
L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).
Le test d’évaluation de la significativité du lien linéaire entre les deux variables est valide, si les résidus :
On vérifie donc en premier l’hypothèse d’indépendance des résidus. On peut représenterl’autocorrélation des résidus à l’aide d’un lag_plot
bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_bw()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
Nous pouvons observer qu’une auto-corrélation significative est présente pour les lags 1,2,3,5,7,10 et 12. Il peut donc y avoir une auto-corrélation pour un lag de valeur 1. Nous pouvons utiliser le test de Durbin-Watson pour le confirmer.
durbinWatsonTest(lin_model_edu)
## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.004
## Alternative hypothesis: rho != 0
Ici, le test nous indique qu’il existe une auto-corrélation significative entre les résidus d’une ligne du tableau de données et ceux de la ligne suivante. Cela signifie qu’a moins que certain residut agissent par groupe, le test rejette le lien linéaire entre les deux variables.
Ensuite on verifie donc que les résidut sont distribués selon une loi Normale de moyenne 0. Cette hypothèse peut s’évaluer graphiquement à l’aide d’un QQplot. Si les résidus sont bien distribués le long de la droite figurant sur le plot, alors l’hypothèse de normalité est acceptée.
qq_edu <- lin_model_edu %>%
ggplot(aes(sample = rstandard(.))) +
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
qq_edu <- ggplotly(qq_edu)
qq_edu
Le test de Shapiro-Wilk peut également être employé pour évaluer la normalité des résidus. L’hypothèse de normalité est rejetée si la p-value est inférieure à 0.05.
shapiro.test(residuals(lin_model_edu))
##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406
Enfin nous testons que les residus sont distribués de façon homogènes. pour cela il faut réaliser un “residuals vs fitted plot”. Les “fitted” correspondent aux réponses prédites par le modèle, pour les valeurs observées de la variable prédicitive.
sc_loc_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) +
geom_point(size = 1.2) +
geom_smooth(method = loess, se = FALSE) +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(0, 1.5),
breaks = seq(0, 1.5, 0.5)
) +
labs(title = "Scale-location",
x = "Fitted values\nlm(prestige ~ education)",
y = "sqrt |Standardized residuals|")
sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu
La courbe rouge, qui est aussi une courbe de régression locale, est globalement plate. Ceci montre que les résidus ont tendance à être répartis de façon homogène tout le long du gradient des valeurs de prestige prédites. Et donc que l’hypothèse d’homogénéité des résidus est acceptée.
Il est également possible d’évaluer cette hypothèse en employant le test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05.
ncvTest(lin_model_edu)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635
Ici, le test ne rejette pas non plus l’hypothèse d’homogénéité.
Il existe aussi la possibilité de directement tester la linéarité à l’aide du plot suivant.
res_v_fit_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), residuals(.))) +
geom_point() +
stat_smooth(method="loess", se = FALSE) +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(-30, 20),
breaks = seq(-30, 20, 10)
) +
labs(title = "Residuals vs. Fitted values",
x = "Fitted values\nlm(prestige ~ education)",
y = "Residuals")
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu
Ceci nous montre que lorsque les réponses prédites par le modèle (fitted values) augmentent, les résidus restent globalement uniformément distribués de part et d’autre de 0. Cela montre, qu’en moyenne, la droite de régression, est bien adaptée aux données, et donc que l’hypothèse de linéarité est acceptable.
summary(lin_model_edu)
##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16
confint(lin_model_edu) %>% kable()
| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | -18.027220 | -3.436744 |
| education | 4.702223 | 6.019533 |
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()
| fit | lwr | upr |
|---|---|---|
| 38.64170 | 36.58966 | 40.69374 |
| 48.61293 | 46.81134 | 50.41452 |
| 67.64405 | 64.52387 | 70.76423 |
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()
| education | income | women | prestige | census | type | residuals |
|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 |
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()
| education | income | women | prestige | census | type | residuals | fitted |
|---|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 | 59.54913 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 | 54.99238 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 | 57.72643 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 | 50.48924 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 | 67.64405 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 | 73.11215 |
vcov(lin_model_edu) %>% kable()
| (Intercept) | education | |
|---|---|---|
| (Intercept) | 13.520978 | -1.1835055 |
| education | -1.183506 | 0.1102162 |
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)
linear_regression <- data %>%
ggplot(aes(y = prestige, x = education)) +
geom_point(size = 1.2) +
geom_smooth(color = "red", method = "lm", fill = "red") +
geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
annotate("text", x = 9, y = 80,
label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3),
" + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(25, 75, 25)
) +
labs(title = "Linear Regression",
x = "Education", y = "Prestige")
linear_regression <- ggplotly(linear_regression)
linear_regression